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Abstract 

We study the statics and dynamics of dark solitons in a cigar-shaped Bose-Einstein condensate 
confined in a double-well potential. Using a mean-field model with a non-cubic nonlinearity, ap- 
propriate to describe the dimensionality crossover regime from one to three dimensional, we obtain 
branches of solutions in the form of single- and multiple-dark soliton states, and study their bifur- 
cations and stability. It is demonstrated that there exist dark soliton states which do not have a 
linear counterpart and we highlight the role of anomalous modes in the excitation spectra. Partic- 
ularly, we show that anomalous mode eigenfrequencies are closely connected to the characteristic 
soliton frequencies as found from the solitons' equations of motion, and how anomalous modes are 
related to the emergence of instabilities. We also analyze in detail the role of the height of the 
barrier in the double well setting, which may lead to instabilities or decouple multiple dark soliton 
states. 
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I. INTRODUCTION 



For many purposes atomic Bose-Einstein condensates (BECs) can be accurately described 
by means of a mean-field theoretical model, the Gross-Pitaevskii (GP) equation l|-|3|. Im- 
portantly the GP equation describes not only the ground-state, but also macroscopic non- 
linear excitations of BECs, such as matter-wave dark solitons, which have been studied 
theoretically (see, e.g., Refs. js, 4| and references therein) and were observed in a series of 



experiments 



5l-llq|. In fact, these localized nonlinear structures have attracted much at 



16 



, [l7|L, while 
. Ad- 



tention as they arise spontaneously upon crossing the BEG phase-transition 
their properties may be used as diagnostic tools for probing properties of BEGs 
ditionally, applications of matter-wave dark solitons have been proposed: the dark soliton 
position can be used to monitor the phase acquired in an atomic matter-wave interferometer 



in the nonlinear regime 19 



As matter-wave dark solitons are known to be more robust in the quasi one-dimensional 
(ID) geometry, the majority of relevant theoretical studies have been performed in the 
framework of the ID GP equation and, particularly, in the so-called Thomas- Fermi (TF)-ID 
regime (see, e.g., Ref. 13|). More specifically, many works have been devoted to the stability 
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- |23| and dynamical properties of dark solitons, such as their oscillations 2lN30| and sound 



emission 
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3l| in the presence of a harmonic trapping potential. In the TF-ID regime, 
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matter-wave dark solitons have also been studied in periodic (optical lattice) potentials 
36], as well as in combinations of harmonic traps and optical lattices j37-4o|. Notice that 
upon properly choosing the harmonic trap and optical lattice parameters, it is possible to 



create a double-well potential, as was demonstrated experimentally in Ref. 



4l|. It may also 



be formed by combining a harmonic trap with a repulsive barrier potential, induced by a 



blue-detuned laser beam 



42| . Double- well potentials have been studied in theory 43N50|. 



while they have played a key role in important experimental observations: these include 
Josephson oscillations and tunneling (for a small number of atoms) or macroscopic quantum 
self-trapping and an asymmetric partition of the atoms between the wells (for sufficiently 
large numbers of atoms) 4l|, as well as nonlinear matter- wave interference leading to the 



formation of matter-wave vortices 
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52] and dark solitons 13l. Il5|. 



So far, there exist only a few studies on matter- wave dark solitons in double-well potentials 
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54| . In particular, Ref. (53| was devoted to the study of nonlinearity-assisted quantum 



tunneling and formation of dark solitons in a matter-wave interferometer, which is reahzed 
by the adiabatic transformation of a double-well potential into a single-well harmonic trap. 
On the other hand, in Ref. 5J|, the stability of the first excited state of a quasi-lD BEG 



trapped in a double- well potential was studied and regimes of (in) stability were found; note 
that in the nonlinear regime, the first excited state of the BEG is nothing but the stationary 



"black" soliton — alias "kink" — solution of t 
it is relevant to mention that both Refs. 53 



le pertinent ID GP equation. At this point. 



54j follow the majority of the theoretical 



works on dark solitons in BEGs, as they have also been performed in the TF-ID regime. 
Nevertheless, this regime has not been practically accessible in real experiments: in fact. 



in most relevant experiments the condensates were usually elongated (a 
three-dimensional (3D) objects, while only two recent experiment s Il3l . 



in the so-called dimensionality crossover regime from ID to 3D 



ias "cigar-shaped") 
were conducted 



55|. The observations of 



;hese experiments were found to be in very good agreement with the theoretical predictions 



13 
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see also Ref. [56|, based on the use of effective ID mean-field models devised in 
Refs. id-y]). 

In the present work we study dark soliton states in cigar-shaped BEGs — being in the 
dimensionality crossover regime from 3D to ID — confined in double- well potentials that are 
formed by a combination of a harmonic trap and an optical lattice. Our analysis relies on the 
use of an effective ID mean-field model, which was first introduced in Ref. 57| and later was 



rederived and analyzed in detail in Refs. (59j . This model, which has the form of a GP-like 
equation with a non-cubic nonlinearity, was also successfully used in Ref. 15| to analyze dark 
soliton statics and dynamics observed in the experiment. Our aim is to systematically study. 



apart from the single stationary dark soliton (see Ref. 54] for a study in the TF-ID regime 



multiple dark soliton states as well; such a study is particularly relevant, as many dark 



solitons can be experimental 
as demonstrated in Refs. 13 



ly created by means of the matter-wave interference method, 



15|. Our study starts from the no n- interacting (linear) limit. 



where we obtain the pertinent excited states of the BEG and, then, using continuation in 
the chemical potential (number of atoms), we find all branches of purely nonlinear solutions, 
including the ground state and single or multiple dark solitons, and their bifurcations. An 
important finding of our analysis is that the optical lattice (which sets the barrier in the 
double-well setting) results in the emergence of nonlinear states that do not have a linear 
counterpart, such as dark soliton states, with solitons located in one well of the double- well 



potential. As concerns the bifurcations of the various branches of solutions, we show that 
particular states emerge or disappear for certain chemical potential thresholds, which are 
determined analytically, in some cases, by using a Galerkin-type approach 5^. For each 
branch, we also study the stability of the pertinent solutions via a Bogoliubov-de Gennes 
(BdG) analysis, which reveals the role of the anomalous (negative energy) modes in the 
excitation spectra as concerns the emergence of oscillatory instabilities of dark solitons: it 
is found that such instabilities occur due to the collision of an anomalous mode with a 
positive energy mode. Furthermore, based on the weakly-interacting limit of the model, we 
study analytically small-amplitude oscillations of the one- and multiple-dark-soliton states. 
Particularly, we obtain equations of motion for the single and multiple solitons from which 
we determine the characteristic soliton frequencies; the latter, are found to be, in many cases 
in good agreement with the eigenfrequencies found in the framework of the BdG analysis. 

The paper is organized as follows. We present our model in section II and provide the 
theoretical framework of our analysis. In Sec. Ill, we study the bifurcations of the branches 
of the solutions of the model and study their stability via the BdG equations. In Sec. IV, 
we provide some analytical estimates for the soliton frequencies based on the dynamics of 
solitons and compare our findings to the ones obtained by the BdG analysis. Sec. V is 
devoted to a study of the soliton dynamics numerically, and shows the manifestation of 
instabilities when they arise. Finally, in Sec. VI, we present our conclusions. 



II. MODEL AND THEORETICAL FRAMEWORK 



A. The effective ID mean-field model and BdG analysis 



We consider a BEG confined in a highly elongated trap, with longitudinal and transverse 
confining frequencies (denoted by and respectively) such that Uz ^ a;_L. In such a 
case, use of the adiabatic approximation, together with a variational approach for determin- 



ing the local transverse chemical potential. 



longitudinal condensate's wave function 
ihdtip ■ 



eads to the following ID GP equation for the 



59| 



ext I, 



+ na;xVl + 4a| 



(1) 



where i^i^Zyt) is normalized to the number of atoms, i.e., = \ip\'^dx, a is the s-wave 
scattering length, m is the atomic mass, and Vext{z) is the longitudinal part of the external 



trapping potential. As demonstrated in Refs. 59|, Eq. ([T]) provides accurate results in 
the dimensionality crossover regime and the TF limit, thus describing the axial dynamics 
of cigar-shaped BECs in a very good approximation to the 3D GP equation. It is worth 
mentioning that in the weakly- interacting limit, 4a|?/'p ^ 1, Eq. ([T]) is reduced to the usual 
ID GP equation with a cubic nonlinearity, characterized byan effective ID coupling constant 
giD = 2ahu± (see, e.g., discussion and references in Ref. |3|). 

Let us now assume that the trapping potential is a superposition of a harmonic trap and 
an optical lattice, characterized by an amplitude Vq and a wavenumber k, namely, 

Ve^t{z) = ^mu^z^ + Vo cos^kz). (2) 

It is clear that a proper choice of the potential parameters (and the number of atoms) may 
readily lead to a trapping potential that has the form of an effective double well potential; 
in such a setting, the height of the barrier can be tuned by the amplitude of the optical 
lattice. In our considerations below, we consider the case of a ^''Rb condensate, confined in 
a harmonic trap with frequencies uj± = lOuz = 27r x 400Hz (these values are relevant to the 



experiments of Refs. 13|, ll5|); furthermore, we will vary the optical lattice strength in the 
range of Vq = to Vq = 1.16 x 10^^^ eV, thereby changing the trapping potential from a 
purely harmonic form to a double-well potential. Finally, we will assume a fixed value of 
the optical lattice wavenumber, namely k = 7r/5.37 fim~^. 

Next, assuming that the density, length, time and energy are measured, respectively, in 



units of a, aj_ = ^JTifmhj\_ (transverse harmonic oscillator length), UJ^ ^ and fux)^^ we express 
Eq. ([T]) in the following dimensionless form. 



idt^ = Hi, + v/l + 4|^|2z^, (3) 

where H = —{l/2)d1 + (1/2)^2^^^ -|- Vocos2(/cz), is the usual single-particle Hamiltonian 
and Vt = uJz/u}± is the normalized harmonic trap strength (note that the normalized optical 
lattice strength Vq is measured in the units of energy hj±). 

Below, we will analyze the stability of the nonlinear modes of Eq. ([3]) by means of the 
Bogoliubov-de Gennes (BdG) equations (see, e.g., Ref. [l|). In particular, once a numerically 
exact — up to a prescribed tolerance — stationary state, ipoiz), is found (e.g., by a Newton- 
Raphson method), we consider small perturbations of this state of the form, 

ij{z, t) = [M^) + {u{z)e-'-' + v*{z)e'-'')] e-'^\ (4) 
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where * denotes complex conjugate. This way, we derive from Eq. ([3]) the following BdG 
equations: 

[H — ji + f]u + gv = uu, (5) 
[H - fi + f]v + gu = -ujv, (6) 

where /i is the chemical potential, and the functions / and g are given by / = + \/l + 4?7,o, 
and g = 2'?/'q/a/1 + Auq (with uq = \ipo\^). Then, solving Eqs. ([5])-([6]), we determine the 
eigenfrequencies u = Ur + iooi and the amplitudes u and v of the normal modes of the 
system. Note that if u is an eigenfrequency of the Bogoliubov spectrum, so are —uj, u* 
and —u* (due to the Hamiltonian nature of the system), and consequently the occurrence 
of a complex eigenfrequency always leads to a dynamic instability. Thus, a linearly stable 
configuration is tantamount to Ui = 0, i.e., all eigenfrequencies being real. An important 
quantity resulting from the BdG analysis is the amount of energy carried by the normal 
mode with eigenfrequency u, namely 

E = [ dz{\uf - \vf)u. (7) 



The sign of this quantity, known as Krein sign [60|, is a topological property of each eigen- 
mode. Importantly, if the normal mode eigenfrequencies with opposite energy (Krein) signs 
become resonant then, in most cases, there a ppe ar complex frequencies in the excitation 



spectrum, i.e., a dynamical instability occurs 



60j. In order to further elaborate on such a 



possibility, it is relevant to note that modes with complex or imaginary frequencies carry 
zero energy (due to the condition {u — u*) J dzdul"^ — |t;p) = 0, which must hold for each 
BdG mode), while anomalous modes have negative energy (see, e.g.. Sec. 5.6 of Ref. [l|). 
The presence of anomalous modes in the excitation spectrum is a direct signature of an 
energetic instability, which is particularly relevant to the case of dark sohtons {2^; such 
excited states of the system are discussed below. 



B. Dark soliton states 



Exact analytical dark soliton solutions of NLS equations with a generalized defocusing 
nonlinearity, such as Eq. (|3]), are not available. In fact, in such cases, dark solitons may only 
be found in an implicit form (via a phase-plane analysis) or in an approximate form (via the 
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small-amplitude approximation) 6l|]. Nevertheless, exact analytical dark soliton solutions 
of Eq. ([3]) can be found in the weakly- interacting limit (4|?/^p <^ 1), and in the absence of the 
external potential: in this case, Eq. is reduced to the completely integrable cubic NLS 
model, which possesses single- and multiple-dark-soliton solutions on top of a background 
with constant density n = nn_= fi (with /i being the chemical potential). A single dark 



soliton solution has the form 



62| 



ip{z, t) = y/no [iu + B tanh(?7)] exp(— Zyut), 



(8) 



where t] = y^B[z — zo{t)], zo{t) is the soliton center, the parameter B = \/\ — v'^ sets 
the soliton depth given by y/noB, while the parameter u sets the soliton velocity, given 
by dzo/dt = ^Jn^u. Note that for z/ = the dark soliton becomes a black soliton (alias a 
stationary kink), with a tt phase jump across its density minimum. Aside from the single- 
dark-soliton, multiple dark soliton solutions of the cubic NLS equation are also available 
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-l64j. In the simplest case of a two-soliton solution, with the two solitons moving with 
opposite velocities, vx = —1^2 = z/, the wave function can be expressed as [6J| (see also 



Ref. 



65|) 



F{z,t) 



exp{—ifit), 



where F 



G{z,t) 

2{no — 2nmin) cosh(T) — 2noi^cosh(Z) + isinh(T), G 



(9) 

2^710 cosh(T) + 



2^ninin cosh(Z'), while Z = 2y/noBz, T = 2A/nmin(no — 'n.min)^, and rimin = no—riQB^ 



is the minimum density (i.e., the density at the center of each soliton). 

Generally, the single dark soliton, as well as all higher-order dark soliton states, can be 
obtained in a stationary form from the non-interacting (linear) limit of Eq. ([3]), corresponding 



to iV ^ 
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67|. In this limit, Eq. ^ is reduced to a linear Schrodinger equation for 



a confined single-particle state, namely the equation of the quantum harmonic oscillator, 
together with the contribution of the optical lattice. However, due to the presence of the 
harmonic trap, the problem is characterized by discrete energy levels and corresponding 
localized eigenmodes. In the weakly-interacting limit, where Eq. ([3]) is reduced to the cubic 
NLS equation, all these eigenmodes exist for the nonlinear problem as well 
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67|, describing 



an analytical continuation of the linear modes to a set_of nonlinear stationary states. Recent 
analysis and numerical results 
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(see also Ref. [15|) suggest that there are no solutions 
of Eq. ([3]) without a linear counterpart, at least in the case of a purely harmonic potential. 
However, as we will show below, the presence of the optical lattice in the model results in 



the occurrence of additional states, for sufficiently large nonlinearity (large atom numbers), 
which do not have a linear counterpart. 

Next, let us discuss the dynamics of dark solitons in the considered setup with an external 
potential in the weakly-interacting limit (4|?/;p ^ 1). Therefore, we factorize the total 
wavefunction ip = UqU into the background density Uq and the soliton wavefunction u. 
We assume that the background density is given by the ground state wavefunction of the 
condensate and that we can use the TF approximation to describe the latter. Then Eq. ([3]) 
may be expressed as a perturbed NLS equation for the soliton wavefunction u, namely, 

idtu + ^d^u- n{\u\'^ -l)u = P{u), (10) 

where the effective perturbation P{u) is given by: 

P(u) = (l-\u\-)Vu+-^-^-g-^. (11) 

Here, V{z) = (l/2)f2^2^ + Vocos'^{kz) is the trapping potential, and all terms of P{u) are 
assumed to be of the same order. Furthermore, assuming that the length scale of the trap is 



much larger than the width of the so 
for dark solitons devised in Ref. 25 



iton, one can apply the adiabatic perturbation theory 



69| to derive the following equation of motion for the 



slowly- varying dark soliton center ZQ{t) 40| : 



df^ 2 dzo ' 
where the effective potential felt by the soliton is given by: 



(12) 



7r2 \ k^^ 



^ 3 y 



cos{2kz). (13) 



At this point we should make the following remarks. First, for the derivation of Eq. ( fT2ll . 
it was assumed that the background density of the condensate can be described via the TF 
approximation; this actually means that Eq. ( IT2|) is valid only for sufficiently large number 
of atoms. At the same time, however, the derivation of Eq. f|T2l) was done in the framework 
of the cubic NLS model, i.e., the weakly-interacting limit of Eq. which is valid for 
sufficiently small number of atoms. Nevertheless, Eq. ( fT2l) may still be relevant to provide 
some estimates. It is known that the soliton oscillation frequency (for a BEC confined in 
a harmonic trap) is up-shifted for large number of atoms due to the dimensionality of the 
system (i.e., the effect of transverse dimensions, which are taken into regard in the derivation 
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of Eq. ([3])) 561] • Particularly, one should expect that in the case of a purely harmonic trap 



the soliton oscillates with a frequency Wds > f^/"\/2j which can be found numerically in 
the absence of the optical lattice (note that Vl/\/2 is the characteristic value of the soliton 



oscillation frequency in a ID harmonic trap of strength Vt in the TF limit |2lN29|). For 
a more detailed discussion see Sec. IIV A[ Moreover, for sufficiently small optical lattice 
wavenumbers (such as the chosen value oi k = 7r/5.37/im~^), and large chemical potentials 
(TF limit) the effective potential of Eq. f|T3|) can be simplified to the following form: 

V,s{z) = culz^ + ^Vo cos{2kz). (14) 

We complete this Section by mentioning that the case of two (or more) spatially well sep- 
arated solitons can also be treated analytically, taking into account the repulsive interaction 
potential between two solitons. In the case of two spatially separated solitons located at zi 
and Zo, on top of a constant background with density no, the interaction potential takes the 



form 



15| 



V-\z,, z.) = . ■ ^. (15) 

smh {y/no{z2 - Zi)) 



III. BIFURCATION AND BDG ANALYSIS 

Let us now proceed by investigating the existence, the stability and possible bifurcations 
of the lowest macroscopically excited states of Eq. ([3]). We fix the parameter values as 
follows: Q = 0.1, k = 7r/5.37 /im~^ and Vq = 1.16 x 10"^^ eV; for such a choice, the four 
lowest eigenvalues of the operator 1 + H (corresponding to the linear problem) are found to 
be: 

uo = 2.138 X IQ-^^eV, ui = 2.142 x IQ-^^eV, U2 = 2.634 x IQ-^^eV, = 2.700 x lO'^^eV, 

(16) 

and the energy of the first excited state is almost degenerate with that of the ground state, 
i.e., ~ uji. 

Fig. [T] shows the number of atoms as a function of the chemical potential for different 
branches of solutions. The insets show the spatial density profiles for /i = 3 x 10~^^ eV, with 
the units of the horizontal and vertical axes being given by fim and 1/fim, respectively. 

Branch (1) corresponds to the states with the largest number of atoms, for fixed chemical 
potential. The respective states are symmetric and have no nodes. Continuation of this 




FIG. 1. (Color online) Number of atoms as a function of the chemical potential for the different 
states. The potential parameters are $7 = 0.1, k = 7r/5.37 fim"^ and Vq = 1.16 x 10^^^ eV. The 
insets show the densities of the different states for /i = 3 x 10~^^ eV. 



symmetric branch to the hnear hmit ends at the eigenvalue cuq, corresponding to the ground 
state of the hnear problem. 

The second branch [branch (2)] starts from the first excited state in the linear limit, 
/i — > cJi as — > 0. The wavefunctions of the states of this branch are antisymmetric and 
have a node at the center of the barrier, while their densities are symmetric. From this 
density symmetric one-soliton branch, two asymmetric one-soliton branches, namely branch 
(3) and its mirror image with respect to the z = axis, bifurcate close to the linear limit, 
at /i ^ 2.144 X 10"^^ eV - see the close-up inset. Let us define a local occupation number, 
i.e., number of atoms in the different wells, as the integral over the density up to the center 
of the barrier (see also Sec. V below). The occupation number in the well with the node 
is then smaller than the occupation number without the node. This can be explained by 
the fact that the state with one node is the first excited state, characterized by a higher 
energy than the one without a node, which is the ground state. Thus, in order to balance 
the chemical potential in both wells one needs a larger interaction energy (i.e., more atoms) 
in the well without a node. Note that the macroscopic quantum self-trapping (MQST) 



state as predicted in Ref 43| is not considered in what follows. The reason for that is that 



the MQST is a running phase state while our bifurcation analysis below focuses on the 
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stationary states of the system. Branches (4) and (5) have no hnear counterparts since at 
fi ~ 2.828 X 10~^^ eV they "coUide" and disappear. The states that belong to these branches 
are asymmetric, exhibiting two nodes. The state with larger number of atoms (for a fixed 
chemical potential) has one node approximately at the barrier and one node in one well. 
The other state has both nodes in one well. Once again, the occupation number in the well 
with two nodes is less than the one of the well with one node, so as to balance the chemical 
potential. 

Branch (6) starts from the second excited state of the linear problem, hence yU — a;2 as 
— 0. Therefore, the density of the respective state is symmetric and there is one node in 
each well. 

Branch (7) starts from the third excited state of the linear problem, and /i — )■ W3 as 

— 0. The states belonging to this branch have three nodes, one located in each well and 
one at the barrier, and are symmetric. Close to the linear limit, at /i ^ 2.795 x 10~^^ eV, 
two asymmetric three-soliton branches, namely branch (8) and its mirror image with respect 
to the z = axis, bifurcate from this state. This state has two nodes in one well and one 
in the other. The occupation number in the well with more nodes is smaller than in the 
well with just one node in order to balance the chemical potential. Notice that there exist 
two more three-soliton branches without linear counterparts, but they only occur at higher 
chemical potentials — where the BEC has occupied four-wells rather than two — so they 
will not be considered here. 

Next, let us study the stability of the states belonging to the above mentioned branches 
by considering the respective excitation spectra, shown in Fig. |2J In this figure, the (blue) 
* symbol denotes a positive Krein sign mode (or a zero one for a vanishing frequency), 
the (green) x symbol a negative Krein sign mode, i.e., a negative energy anomalous mode, 
and the (red) (+) symbol a vanishing Krein sign, i.e., an eigenmode associated with com- 
plex/imaginary eigenfrequency. Generally speaking, there are two different kinds of insta- 
bility, corresponding to the cases of either a purely imaginary eigenfrequency, or a genuinely 
complex eigenfrequency. The latter case gives rise to the so-called oscillatory instability, 
stemming from the collision of a negative energy mode with one of positive energy. 

The first panel of Fig. |2] shows the BdG spectra along the first branch. The imaginary 
part is zero for every value of the chemical potential, a fact reflecting the stability of this 
state. The absence of negative energy modes is expected, as this state is actually the ground 
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H(10"'^eV) n(10"'^eV) |i(10"'^eV) |i(10"'^eV) 



FIG. 2. (Color online) BdG Spectrum of the states shown in Fig. [TJ First panel shows the 
spectrum of the first branch, second panel the spectrum of the second branch and so forth. The 
(blue) symbol denotes an eigenmode of positive Krein sign (or zero Krein sign for a vanishing 
eigenfrequency) , the (green) x symbol an eigenmode of negative Krein sign, and the (red) + symbol 
an eigenmode associated to a vanishing Krein sign with complex/imaginary eigenfrequency. 

state of the system. Since the BdG spectrum refers to all excitations of the state, one can 
attribute in the linear limit to each mode an excited state of the linear problem. The energy 
gap between the ground state and the first excited state corresponds to the energy of the 
first mode. Since the frequency of the so-called bosonic Josephson junction, namely the 
oscillation frequency characterizing the transfer of atoms between the two wells, is equal to 
this frequency in the linear limit, this mode is connected to the oscillation of atoms between 
the two wells. We investigate this fact in more detail below (see Sec. |V]) by direct simulations 
in the framework of Eq. ([T]). The frequency value of the second nonzero BdG mode of the 
ground state spectrum is correspondingly equal to the energy gap between the ground state 
and the second excited state in the linear limit, and so forth. We define the difference in 
the way that these frequency differences are positive. 
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The second panel shows the BdG spectra along the second branch, solutions of which 
represent states with a density minimum at the center of the barrier. In the linear limit 
one can assign to each mode an energy difference of the linear problem similar as in the 
BdG spectrum of the ground state. However in this case the frequency difference of the 
first excited state to the ground state is negative (according to our definition) leading to a 
negative energy mode. As one can see in the inset, the lowest mode has negative energy close 
to the linear limit. However, increasing the chemical potential the eigenfrequency of this 
mode is moving rapidly towards the origin of the spectrum and, at /i = 2.1437 x 10^^^ eV, 
becomes and remains purely imaginary, thus carrying zero energy (red + symbols). This 
happens exactly at the point where the third branch bifurcates as one can see in the inset 
of Fig. [U For the above mentioned potential parameters, this bifurcation happens in the 
weakly-interacting regime, where Eq. (3) reduces to the cubic NLS equation. Thus, one can 



apply the Galerkin-type approach of Ref. 50| and find that the critical value of the norm 



at which the bifurcation occurs is given by 

where wq,! are the first two lowest eigenvalues of the operator 1 + H, Ai = J ipfdz, B = 
j ipQipfdz, while denote the ground state and the first excited state of the operator 
1 + i^^, respectively. One can also determine the critical chemical potential at which the 



symmetry breaking is expected to occur, namely ficr = 1 + + ^iBN^r 50| leading for (these 
parameter values) to ficr = 2.1436 x 10"^^ eV, which is in very good agreement with our 
numerical findings. 

The third panel shows the BdG spectra along the third branch which bifurcates at = 
2.1437 X 10"^^ eV from the second one. All the eigenfrequencies of the arising asymmetric 
state are real close to the linear limit, thus it is concluded that this state is linearly stable 
and that the bifurcation is a supercritical pitchfork. This state has one negative energy 
mode as well, since it is a one-soliton state. The eigenfrequency of the negative energy mode 
increases with increasing chemical potential, passes through the mode corresponding to the 
second excited state of the linear limiting case, and finally collides at some critical value 
of the chemical potential, /i = 2.76 x 10~^^ eV, with the mode corresponding to the third 
excited state; this results in the generation of a quartet of unstable eigenfrequencies and, at 
this point, the state becomes oscillatory unstable. However, the magnitude of the instability 
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is much smaller than in the previous case of the purely imaginary eigenfrequency. 

The fourth and fifth panel show the spectra along the fourth and fifth branch, the solutions 
of which correspond to asymmetric two-soliton states having no linear counter-part. The 
previous one has got one imaginary and one negative energy mode and the latter one has 
two negative energy modes. These modes result from the presence of two dark solitons. The 
state considered in the fourth panel has two dark solitons, one soliton located approximately 
at the barrier and one soliton in one well. The eigenfrequency corresponding to the former 
one is purely imaginary and looks similar to the mode of the symmetric one soliton state of 
the second panel. The eigenfrequency of the other anomalous mode decreases with increasing 
chemical potential and collides with a mode with positive energy, thus generating a quartet 
of complex eigenfrequencies. The magnitude of the imaginary part is again much smaller 
than the magnitude of the purely imaginary eigenfrequency. The generation of dynamic 
instabilities is illustrated e.g. in panel five. The collision of one of the anomalous modes 
with a mode of positive Krein sign could lead to the emergence of a quartet of complex 
eigenfrequencies, resulting in the concurrent presence of an apparent level crossing in the real 
part of the eigenfrequency along with a nonzero imaginary part of the relevant eigenfrequency 
mode. 

At /i ~ 2.828 X 10"^^ eV, the states of the fourth and fifth branch collide and disappear 
as can be seen in Fig. [T] Just before the collision (for slightly larger /i), the state of the 
fourth branch is unstable due to the presence of a mode with imaginary eigenfrequency in 
the BdG spectrum, while the state of the fifth branch is linearly stable since all the BdG 
eigenfrequencies have only real parts. Thus, it can be concluded that, at this critical point, 
a saddle-node bifurcation occurs, readily destroying these two states. 

The sixth panel shows the BdG spectra along the sixth branch, the solutions of which 
correspond to symmetric two-soliton states. This branch has a linear counterpart, and we 
can assign to each mode in the linear limit an energy difference of the linear problem. The 
negative energy modes occur for negative energy differences, e.g., at the energy difference to 
the ground state and the first excited state. The anomalous modes here are close to being 
degenerate. Physically speaking, this reflects the fact that the dark solitons are essentially 
decoupled at the linear limit due to the presence of the optical lattice. The negative energy 
modes are not the modes with the lowest eigenfrequency. The energy gaps between the 
second excited state and the third and fourth excited state, respectively, are smaller than 
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the gaps between the ground state and the first excited state. Therefore, these modes with 
positive energy occur in the linear hmit at lower frequencies. One of the negative energy 
modes collides at some critical value of the chemical potential with the mode stemming 
from the fourth excited state and forms a quartet of complex frequencies, thus leading to an 
oscillatory instability. 

The seventh panel shows the BdG spectra along the seventh branch, the solutions of 
which corresponds to the symmetric three-soliton state exhibiting a linear counterpart. The 
negative energy modes occur at the energy gaps to the ground state, the first and the second 
excited states. The eigenfrequency of the latter becomes, for increasing fi, purely imaginary 
at the point where the eighth branch bifurcates. The other two anomalous modes are 
almost degenerate, reflecting again the fact that the solitons in the different wells are almost 
decoupled. One of the negative energy modes collides with a mode with positive energy 
and forms a quartet of complex frequencies. However, the magnitude of the imaginary part 
is much smaller than the magnitude of the purely imaginary mode. The behavior of the 
purely imaginary mode is similar to the behavior of the mode in panel two (describing a 
single soliton at the center of the barrier). The behavior of the other two modes is similar 
to the modes in panel six describing one soliton in each well. 

Finally, the eighth panel shows the BdG spectra along the eighth branch, the solutions of 
which corresponds to the asymmetric three-soliton state. This state, which has two solitons 
in one well and one soliton in the other well, emerges at ~ 2.74 x 10~^^ eV, where the 
seventh state (i.e., the one corresponding to the nonlinear continuation of the third excited 
state), becomes unstable. The critical value of the number of atoms for which the bifurcation 
occurs can be predicted in the same way as for the bifurcation of the one soliton branch; 
the result is: 

(3) ^ ^3-C^2 , . 

" 3D - Ci ' ^ ^ 

where 072,3 are the third and fourth lowest eigenvalues of operator 1 + H, Ci = J ip^dz, 
B = J ipli'ldz, while V^2,3 denote the second and third excited state of the operator 1 + if, 
respectively. The corresponding critical chemical potential is then given by jjL^cr = 1 + ^2 + 
3BNcr^ leading (for these parameters) to /ii?'* = 2.7420 x 10~^^ eV, which is in very good 
agreement with the numerical critical point given above. 

Right after the bifurcation, the asymmetric three soliton states are stable and, thus, we 
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can conclude that the bifurcation is of the supercritical pitchfork type. Two modes look 
similar to the negative energy modes in panel five describing a two-soliton state in one well 
and the third mode looks similar to a mode of one-soliton in a harmonic trap. 

From the above analysis, we can readily derive some general conclusions concerning dark 
soliton states in a double well potential. The sum of negative energy modes and imaginary 
eigenfrequency modes is equal to the number of nodes in the wavefunction profile which, for 
large values of the chemical potential, represent dark solitons. In the linear limit, the BdG 
eigenfrequencies of a state correspond to the energy differences between that and other states 
of the linear system. For a negative energy difference the corresponding mode has a negative 
energy. Dark solitons located at the center of the barrier are known to be linearly unstable 
for sufficiently high chemical potential) associated with a purely imaginary eigenfrequency 
Thus, symmetric states with an odd number of nodes always have a mode which. 



above a critical value of the chemical potential, becomes purely imaginary. At this critical 
point additional stable asymmetric dark-soliton states emerge through supercritical pitchfork 
bifurcations. 



IV. STATICS VS. DYNAMICS OF DARK SOLITON STATES 
A. The one-soliton state 

Having investigated in detail the statics of matter- wave dark solitons in double-well poten- 
tials, we now proceed to compare these results to the soliton dynamics, using the theoretical 
background exposed in Sec. Illll B. Particularly, considering the case of a single dark soliton, 
we can readily obtain from Eqs. ( [T2|) and (fT4|) the following approximate equation of motion 
for the dark soliton center: 

^ = -ulzo + sin{2kz^). (19) 

It is straightforward to find the fixed points Zo^Po = dzo/dt of the above dynamical system, 
which are given by 

Po = 0, Zo = sm{2kzo). (20) 

It is clear that the fixed points correspond to intersections of the straight line fi{zo) = Zo 
and the sinusoidal curve f2{zo) = sm{2kZo). The fixed point at (-^cPo) = (0,0) exists 

ds 
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for all choices of Vq, k). This equilibrium corresponds to the symmetric one-soliton state, 
namely a stationary black soliton located at the center of the potential. In the absence of 
the optical lattice (i.e., for Vq = 0) this is the only existing fixed point, i.e., in the case of 
a purely harmonic trap, only a symmetric one-soliton state can exist, which originates from 
the first excited state of the linear problem. On the other hand, when the optical lattice is 
present (i.e., for Vq 7^ 0), this situation changes and more fixed points arise when the slope of 
the curve f2{zo) at the origin becomes larger than the one of the straight line fi{zo), namely 
when /c^Vo/cjjg > 1, two new fixed points appear through a supercritical pitchfork, with 
the newly emerging states corresponding to the bifurcating asymmetric one-soliton states. 
Therefore, the bifurcating state appears only for a sufficiently strong optical lattice, such 
that Vq > Vo,cr = ^dsl^"^- Note that for even larger values of the parameter /c^Vo/w^g more 
fixed points appear, corresponding to solitons located in further (i.e., more remote from 
the trap center) wells of the optical lattice; for this reason, these states are not considered 
herein. 

Let us now investigate the spectral stability of the fixed points employing the linearization 
technique (see, e.g., Ref. Particularly, we first find the corresponding Jacobian matrix 

(evaluated at the fixed points), which is given by: 



When evaluated at the trivial fixed point (^07^0) = (0,0), the corresponding eigenvalue 
problem for the Jacobian, det(J — AI) = 0, leads to eigenvalues = —oj^^ + k'^Vo. It is clear 
that in the case of /c^Vo/w^g < 1, the eigenvalues are imaginary (note that, accordingly, the 
eigenfrequencies are real since u = iX) and, thus, the fixed point (0, 0) is a center. In this 
case, if the dark soliton is weakly displaced from the center of the trap, it performs harmonic 
oscillations around the center with an oscillation frequency given by Uosc = \/^ds ~ ^^^o- 
On the other hand, if fc^-^ > 1 then the eigenvalues become real (and, accordingly, the 
eigenfrequencies are imaginary), hence (0,0) becomes a saddle point, while the respective 
imaginary eigenfrequency is given by Wosc = i\/\^ds ~^ For the newly bifurcating 

(stable, since Vq > u^^/k"^) fixed points one can obtain the following eigenfrequency 





(21) 
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From the above analysis it is clear that the supercritical pitchfork occurs at 

Vo,cr- = -j^- (22) 

Based on the above, one can make the following estimates and comparisons to numerical 
results: in the TF limit, and in the weakly-interacting case (i.e., in the framework of the 
ID cubic GP equation), the dark soliton oscillation frequency is predicted to be Uz/\/2 = 
28.28Hz. However, as mentioned in Sec. IIIIi B. a larger value of this frequency is expected in 
the framework of Eq. ([T]), due to the effect of the dimensionality of the system. Employing 
a BdG analysis in the case of a pure harmonic trap, we find that the anomalous mode 
eigenfrequency (i.e., the soliton oscillation frequency) decreases with increasing chemical 
potential and asymptotically reaches the value Uds = 30.36Hz (see first panel of Fig. H] 
below). Using this value (and for k = 7r/5.37/im~^), one obtains Vo^cr = 9.6 x 10"^"^ eV. 




0.5 1 0.5 1 

Vo(10'''eV) Vo(10-^=^eV) 



FIG. 3. (Color online) Top panel: Fixed points (lines) of the dynamical system (I19p and position 
of the soliton (open circles) calculated by the Eq. ([T]) with /i = 3.8 x 10~^^ eV. Bottom left (right) 
panel: BdG spectra of the symmetric (asymmetric) one-soliton state for = 3.8 x 10~^^ eV as a 
function of Vq. 

In the top panel of Fig. [3l we show the fixed points of the dynamical system (|T9|l as a 
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function of the optical lattice strength Vq. The dash-dotted (green) line denotes the non-zero 
fixed point calculated from Eq. fl20|) with Uds = ^z/ V^, while the solid (black) line denotes 
the non-zero fixed point calculated from Eq. ( 1201) with Uds = 30.36 Hz. Furthermore, the 
circles denote the position of the node of the nonlinear stationary state — i.e., the position of 
the dark soliton — obtained by numerical integration of Eq. ([T]) with /i = 3.8 x 10~^^ eV (in 
the TP limit). Clearly, at Vo^cr — 9.6x10"^^ eV, the asymmetric one-soliton branch bifurcates 
through a supercritical pitchfork bifurcation; notice the excellent agreement between the two 
approaches. 

The bottom left and right panels of Fig. [3] show the BdG spectra of the symmetric 
and asymmetric one soliton state, respectively, as a function of the optical lattice strength 
(for a fixed chemical potential, fi = 3.8 x 10~^^ eV). For sufficiently small lattice strength, 
Vq < Vo,cr; the Symmetric state is stable, while when the lattice strength Vq is increased, the 
magnitude of the instability grows and the prediction of Eq. (|2T1) deviates from the results 
obtained from the integration of Eq. ([3]). To explain this discrepancy, it is relevant to recall 
that the result of the perturbation theory of Sec. [TllB is valid only if all perturbation terms 
are small and of the same order. Nevertheless, for the symmetric one soliton state — with 
the soliton located at zero — the magnitude of the perturbation term corresponding to 
the harmonic trap (optical lattice) has its smallest (largest) value; thus, for large Vq, these 
perturbations cannot be of the same order. As concerns the asymmetric state, the prediction 
of Eq. (I2TI) is in fairly good agreement to the results obtained from the numerical integration 
of Eq. ([3]). However, Eq. ( 12T1) cannot predict oscillatory instabilities (associated with the 
existence of complex eigenfrequencies) , which occur when the anomalous mode collides with 
a mode of the background condensate. This can readily be understood by the fact that the 
derivation of the equation of motion of the dark soliton relies on the assumption that the 
soliton is decoupled from the background (see, e.g., Refs. 25 1 and for details), which is 
not the case for oscillatory instabilities. 

Next, let us study the continuation of both symmetric and asymmetric one-soliton states 
with increasing chemical potential for various values of the optical lattice strength. 

First, we consider the symmetric one-soliton state and in Fig. |4]we show the BdG spectra 
along the corresponding branch. The first panel shows the case of a pure harmonic trap, i.e., 
Vq = 0, and the dashed line is the prediction from Eq. ( 12T|) using tUds = (^z/V^- It is clear 
that the eigenfrequencies are real, indicating the stability of this state for all values of the 
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FIG. 4. (Color online) BdG spectra for four different values of the optical lattice strength for 
the symmetric dark soliton state (with the soliton located at the trap center): from left to right, 
Vo = 0, Vb = 8.3 X 10"^^ eV, Vq = 1.65 x IQ-^^ gy ^nd Vq = 4.96 x lO"^^ eV. The solid lines 
correspond to the predictions of Eq. (|2ip . Symbol * (blue) denotes a positive (or zero for a 
vanishing eigenfrequency) , symbol x (green) a negative and symbol + (red) a vanishing Krein sign 
mode (for non-zero eigenfrequencies). 

chemical potential. The figure shows the constant eigenfrequency of the dipolar mode [see 
straight (blue) line] at the value of the trap frequency ujz and the anomalous mode bifurcating 
from the dipolar mode in the linear limit. As mentioned above, the eigenfrequency of the 
anomalous mode reaches asymptotically the value uj^s = 30.36 Hz. 

The second panel of Fig. |4] shows the BdG spectrum for Vq = 8.3 x 10"^^ eV. All eigen- 
frequencies are real, so the state is still stable. The dashed line is the prediction by Eq. f l2T|) 
using uJds = uJz/V^, while the solid line is the prediction with lo^s = 30.36 Hz. The dipole 
mode remains approximately constant for such a small optical lattice strength. However, 
the degeneracy of the anomalous mode and the dipole mode in the linear limit disappears. 

The third panel shows the spectra for Vq = 1.65 x 10"^^ eV. In this case, the anomalous 
mode becomes purely imaginary for increasing chemical potential, reflecting the fact that 
the state becomes unstable. Notice that the numerical result and the analytical prediction, 
^ODE = ^v^l^^ds ~ ^^K)!; are in a good agreement for large values of the chemical potential, 
where the TF limit is reached and Eq. f|T9|) is valid. Since Vq > Vo,cr; this state is expected 
to be unstable in the TF limit. 

The fourth panel shows the spectra for Vq = 4.96 x 10^^^ eV. In this case, the anomalous 
mode becomes purely imaginary already at small chemical potentials. Furthermore, the 
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dipole mode eigenfrequency is now significantly modified, upon the reported interval of 
changes of the chemical potential. The analytical prediction agrees again well with the 
numerical BdG result in the TF limit. 
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FIG. 5. (Color online) BdG spectra for Vq = 1.65 x lO'^^ eV (left panel) and Vq = 8.3 x 10"^^ eV 
(right panel) for the asymmetric dark soliton state, bifurcating from the symmetric one. The 
notation of eigenstates is the same as in Fig. HI 



Next, we consider the continuation with increasing chemical potential of the asymmetric 
states bifurcating from the symmetric one when the latter becomes unstable. Fig. [5] shows 
the BdG spectra along the asymmetric one-soliton state for Vq = 1.65 x 10~^^ eV and Vq = 
8.3 X lO"^'^ eV. The analytical predictions of Eq. (!2T|) with u^s = oJz/^/^ and = 30.36Hz, 
respectively plotted with dashed and solid lines are in good agreement with the respective 
anomalous mode eigenfrequencies. In the case of Vq = 8.3 x 10"^'^ eV, the difference between 
the result of Eq. fl^ using Uds = i^z/v^ and Uds = 30.36 Hz is small, showing that the 
dynamics is mostly governed by the potential term stemming from the optical lattice. 



B. Two dark-soliton states 

We now consider the case of the two-soliton state, for which we need to incorporate the 
interaction potential in order to describe their dynamics. Particularly, according to the 
theoretical framework of Sec. HHB, the effective potential describing the dynamics of two 
solitons located at zi and Z2 has the following form: 

Vesizu Z2) = oolizl + zl) + ^(cos(2A;2i) + cos{2kz2)) + . '^t TT- (23) 

2 smh [y/rMy[z2 - Zi)) 
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Using d?Zi/dt^ = — (1/2) dVcs{zi, Z2) /dzi (with i = 1,2) as per Eq. ( IT2|) . it is straightforward 
to derive the following system of equations of motion, 

d'^Zi 1 /2 

-j^ = -Wds^i + 2^Vo sin(2fczi) - 2nl''^ coth.{y/n^{z2 - Zi))csch'^ {y/n^{z2 - Zi)), (24) 
^2^ 1 

= -ujIsZ2 + 2^^o sin(2A;z2) + Srig^^ coth(A/ri^(z2 - zi))csch^(v/ri^(2;2 - ^1)), (25) 
which, accordingly, leads to the following system for the fixed points {zi^\ z^^): 



= ^ sin(2fc4^)) - \- coth(y^(42) - 4^)))csch^(y^(.f - 4^))), (26) 

^'^ds ^ds 



ZA- = ^ sin(2fcz(2)) + ^ coth(^(4^) - 4^)))csch^( V^(4^^ - 4^0). (27) 

^'^ds ^ds 



Thus, the fixed points for the two soliton states depend on the density of the background 
cloud and thereby on the chemical potential. Due to the presence of the external potentials, 
the density of the background is inhomogeneous. In our analysis, we replace the inhomo- 
geneous density, with the density of the background at 2 = using the TF approximation; 
namely, no{z = 0) = ((/i - Vof - l)/4. 

Considering small deviations from the fixed points, i.e., Zi = Zo^ + eriiexp{iuit) (with 
i = 1,2), leads to the eigenvalue problem for the normal modes: 

2 ^ lA[z^^^) + B[z^\zf^) -B[z^^\zf^) \ 

\ -B{z^^\z^y) Aizi'^) + Bizi'\zi'^))'' ^ ^ 

with 

A{x) = cul^ - Vok"^ cos(2A;x) (29) 
B{x, y) = Anl coth?{^/n^{x - y))cf^ch^ {^/n^{x - y)) + 2noCsch''(v^(a; - y)). (30) 

Diagonalizing the above matrix, we obtain the eigenfrequencies of the two anomalous modes 
of the two-soliton states. Notice that since different two soliton states lead to different fixed 
points, the oscillation frequencies for the two solitons will differ as well. Negative eigenvalues 
of this 2x2 system lead to imaginary frequencies describing unstable states; however, as 
mentioned in the previous section, one cannot describe oscillatory instabilities. 

In the top panel of Fig. [6l we show the fixed points of the system (l26l) -(l271) and the zero 
points of the solutions with two nodes of Eq. ([1]) corresponding to the positions of the two 
solitons as a function of Vq (the chemical potential is fixed: = 3.6 x 10^^'^ eV). For small Vq 
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FIG. 6. (Color online) Top panels: Fixed points (lines) of the system (j26p -(|27p and position of the 
solitons (circles for the symmetric two-soliton state, x for the two-soliton branch with one soliton 
near the center and one in one well, and stars for the state with both dark solitons in the same 
well) calculated by the numerical integration of Eq. ([T|) for fj. = 3.8 x 10"^^ eV. Bottom panels: 
BdG spectra of the two-soliton states for fi = 3.8 x 10"^^ eV as a function of Vq; left corresponds 
to the symmetric two-soliton state, middle to the always unstable state with one soliton near the 
center and one in one well, while the right corresponds to the state with both dark solitons in the 
same well. Notice also the agreement of the theoretically predicted modes from Eq. (j28p with the 
numerical ones. 

there exists only the symmetric solution with two nodes. However, for Vo,cr = 3.4 x 10~^^ eV, 
two asymmetric states appear stemming from a saddle-node bifurcation. For all three states, 
the agreement between the fixed points obtained from Eqs. fl26l) - fl27j) and the position of the 
solitons is excellent. 

The bottom panels of Fig. [6] show the BdG spectra for the symmetric and the two 
asymmetric two soliton states as a function of Vq for fixed fj, = 3.8 x 10"^^ eV. The solid 
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lines correspond to the normal modes of Eqs. (12^ - (125|) around the corresponding fixed 
points. For a small optical lattice strength Vq, the two anomalous modes of the symmetric 
states have a different magnitude, which indicates that the solitons are coupled. However, for 
increasing Vq, this difference decreases and, thus, the coupling between the soliton becomes 
weaker. Notice that for an optical lattice strength Vq ~ 3.6 x 10"^^ the two modes become 
of the same order, a fact that indicates that the two solitons become eventually decoupled. 
The prediction of Eqs. f lM|) - fl25|) agrees very well with the numerical results obtained via the 
BdG analysis for the anomalous modes. The second panel shows the BdG spectrum of the 
asymmetric state with one soliton located approximately at the center of the barrier. The 
latter, leads to a purely imaginary mode with an increasing magnitude with increasing Vq. 
The increase of the magnitude of this mode results from the increase of the height of the 
barrier. For large Vq the prediction of Eqs. flM|) - fl25|l deviates from the result obtained in the 
framework of Eq. ([1]): the contributions of the perturbation terms are not of the same order 
and, thus, perturbation theory fails in this case. The third panel shows the BdG spectrum 
for the second asymmetric state with two solitons in one well. The difference between the 
modes is always large, a fact indicating the strong coupling between the two solitons. With 
increasing lattice strength Vq, the eigenfrequencies of the anomalous modes increases too. 
The predictions of Eqs. (124|) -( 125|) agree once again well with the BdG results; nevertheless, 
as in previous cases, Eqs. f l24|) -f l25|) cannot predict oscillatory instabilities. 
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FIG. 7. (Color online) BdG spectra for different values of the optical lattice strength; from left to 
right, Vq = 0,Vq = 8.3 x 10^^^ eV, Vq = 1.65 x IQ-^^ gy and Vq = 4.96 x IQ-^^ ^y. The solid lines 
are the predictions of Eqs. ()24p -(j25 p with a numerically determined oscillation frequency (wds)- 

Next, let us consider a continuation with respect to the chemical potential. Fig. [7] shows 



24 



the BdG spectra, and the predicted oscillation frequencies, for the symmetric two soliton 
state for different values of the optical lattice strength: Vq = 0, Vq = 8.3 x 10~^^ eV, 
Vq = 1.65 X 10"^^ eV and Vq = 4.96 x 10"^'^ eV. In the case of a purely harmonic potential, 
Vq = 0, one observes the constant eigenfrequency of the dipolar mode at u^, as well as the 
eigenfrequency of one anomalous mode bifurcating from the dipolar mode in the linear limit. 
The eigenfrequency at 2ujz is imaginary due to a degeneracy of a negative and a positive 
Krein sign mode. The former one corresponds to the second anomalous mode and the latter 
one to the quadrupole mode. This degeneracy is lifted for increasing chemical potential. 
The predictions of Eqs. f l21|) -f l2^ agree very well with the BdG results for the anomalous 
modes in the case of large chemical potentials. 

For Vq = 8.3 x 10"^^ eV the degeneracy between the anomalous modes and the dipolar and 
the quadrupole mode is lifted in the linear limit. Thus, there exists an energy gap between 
the anomalous mode and the dipolar mode and, more importantly, an energy gap between 
the second anomalous mode and the quadrupole mode. A particularly useful conclusion 
stemming from the above findings is that the two soliton state can be stabilized in a harmonic 
trap by adding a small optical lattice. 

For Vq = 1.65 x 10~^^ eV the system is stable as well. The anomalous modes decrease 
with increasing chemical potential and reach similar values for a large chemical potential. 
This denotes that the two solitons are decoupled for a large chemical potential. This can 
be understood at least in the framework of the effective potential fl23|) : the range of the 
interaction potential scales with the density uq of the background atom cloud and thereby 
with the chemical potential as well. Thus, for a large chemical potential, i.e., a large density, 
the interaction between solitons is strong but of short range too. 

For Vq = 4.96 x 10"^^ eV the anomalous modes have similar values even in the linear 
limit. Thus, the two solitons are already almost decoupled in the linear limit due to their 
large separation. The two anomalous modes collide with the dipolar and quadrupole modes 
for increasing chemical potential, leading to oscillatory instabilities, and bifurcate from these 
modes again for a larger chemical potential. Although the predictions of Eqs. agree 
qualitatively with the BdG results, there is no quantitative agreement due to the crossings 
and collisions of the soliton modes with the modes of the condensate. 
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V. DARK SOLITON DYNAMICS: NUMERICAL RESULTS 



So far, we have examined the existence, stabihty and bifurcations of the branches of 
dark-sohtonic states in the double-well potential setting. We now proceed to investigate the 
dynamics of these states by direct numerical integration of Eq. ([1]), using as initial condi- 
tions stationary soliton states, perturbed along the direction of eigenvectors associated with 
particular eigenfrequencies (especially ones associated with instabilities). To better explain 
the above, we should recall that anomalous modes are connected to the positions and oscil- 
lation frequencies of the solitons; therefore, perturbations of stationary dark soliton states 
along the directions of the eigenvectors of the anomalous modes result in the displacement 
of the positions of the solitons, thus giving rise to interesting soliton dynamics. 

Before proceeding with the dynamics of the (excited) soliton states, let us make a remark 
concerning the ground state of the system. As mentioned above, the value of the lowest BdG 
mode of the ground state is given, in the linear limit, by the energy difference between the 
ground state and the first excited state. However, this energy gap determines the oscillation 
frequency of atoms between the wells of the double-well potential, which suggests that 
the lowest BdG mode is connected to an oscillation between the wells. Let us define the 
atom numbers Nj^ and N^i in the left and right well, respectively, as A^l = {ipl'^dz and 
= /o^°° \4'\'^dz (the total number of atoms is = A'l + A^r). In Fig[8l we show the time 
evolution of A^l (solid line) and (dashed line) for the ground state, when perturbed by the 
eigenvector corresponding to the lowest BdG mode (parameter values are /i = 3 x 10~^^ eV 
and Vo = 1.12 x 10"^^ eV). One can clearly observe the oscillation of the atom numbers in the 



different wells, with a characteristic frequency known as Josephson oscillation frequency 
determined by the magnitude of the respective BdG eigenmode. 

Let us now proceed with the one-soliton states. Fig. [9] shows the time evolution of the 
symmetric one soliton state perturbed by the eigenvector corresponding to the unstable 
mode, for yU = 3 x 10~^^ eV and Vq = 4.96 x 10"^'^ eV. The instability manifests itself as a 
drift of the soliton from the trap center to the rims of the background atom cloud, where 
the soliton is reflected back and forth. One clearly observes that the soliton is located at 
an unstable position, since a small perturbation leads to a large amplitude oscillation of the 
soliton. On the other hand. Fig. [10] shows the time evolution of the asymmetric one soliton 
state perturbed by the eigenvector corresponding to the unstable mode for = 3 x 10~^^ 
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FIG. 8. (Color online) Time evolution of the number of atoms A^l (solid line) and (dashed 
line) in the left and right well of the double-well potential, respectively, for the ground state of 
the system perturbed by the lowest BdG mode. Parameter values are /i = 3 x 10"^^ eV and 
Vo = 1.12 X 10-12 eV. 
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FIG. 9. (Color online) Contour plot showing the space-time evolution of the density of the con- 
densate, which carries the symmetric one-soliton state; dark corresponds to minimal and white 
to maximal density. The soliton is unstable and departs from the trap center after a short time 
period, performing subsequently large amplitude oscillations. Parameter values are // = 3 x IQ-^^ 
eV and Vb = 4.96 x IQ-^^. 

eV for Vo = 8.27 x 10^^^ eV. The corresponding state is subject to an oscillatory instability, 
which results in a growing amplitude of the soliton oscillation. However, in comparison to 
the previous case, the magnitude of the instability is much smaller and, thus, the increase 
of the amplitude is smaller (for the same time scale) as well. 

Next, we consider the multi-soliton states and start with the symmetric two-soliton state 
which, for = 3.3 x IQ-^^ eV and Vq = 8.27 x IQ-^'* eV, is found to be linearly stable. 
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FIG. 10. (Color online) Same as in Fig. [21 but for the asymmetric one-soliton state, for = 3x 10~^^ 
eV for Vn = 8.27 x IQ-^^ ^y. 



Numerical simulations, using as initial conditions this state perturbed by the eigenvectors 
corresponding to both anomalous modes, reveal that the two solitons perform oscillations 
with a constant amplitude and different frequencies. Specifically, the mode with smaller 
amplitude performs an in-phase oscillation (with the two solitons moving towards the same 
direction), whereas the mode with larger amplitude performs an out-of-phase oscillation 
(with the two solitons moving towards opposite directions) with a larger oscillation frequency. 
The existence of two different oscillation frequencies indicates the coupling of the dark 
solitons. This situation changes for a large optical lattice strength (e.g., Vq = 1.65 x 10~^^ 
eV): in this case, the two dark solitons are decoupled, since the numerical simulations show 
that both the in-phase and out-of phase oscillations are characterized by almost the same 
frequency. 

Fig.[TT]shows the time evolution of the symmetric two-soliton state for /i = 3.1 x 10"^^ eV 
and Vo = 4.96 x 10~^^ eV perturbed by the eigenvector corresponding to the anomalous mode 
associated with in-phase (left panel) and by the eigenvector corresponding to the imaginary 
mode out-of-phase (right panel) oscillations of the dark solitons. For these parameters, the 
BdG analysis predicts that the former mode is stable while the latter is unstable. This 
is confirmed by the direct simulations: it is clearly observed that the amplitude of the 
out-of-phase oscillation increases whereas the amplitude of the in-phase oscillation remains 
constant. 

Finally, in Fig. [12] we show the time evolution of the symmetric two soliton state for 
= 2.8 X 10^^^ eV and Vq = 4.96 x 10"^'^ perturbed by the eigenvector corresponding 
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FIG. 11. (Color online) Same as in Fig. [9l but for the symmetric two-soliton state with in-phase 
(left) anomalous mode and out-of-phase (right) imaginary mode perturbations, for /x = 3.1 x 10~^^ 
eV and Vq = 4.96 x 10"^^ eV. The former anomalous eigenmode is stable, while the latter is 
unstable (as concluded also from the BdG analysis). 
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FIG. 12. (Color online) Same as in Fig. \TT\ but for the symmetric two-soliton state, for /i = 
2.8 X 10~^^ eV and Vq = 4.96 x 10^^^ eV. Notice the reversal of stability of the two anomalous 
modes with respect to Fig. [TTJ 



to the anomalous mode with in-phase (left panel) and by the eigenvector corresponding 
to the imaginary mode out-of-phase (right panel) oscillations. For these parameters, the 
simulations show that the out-of-phase oscillation of the solitons is stable and the in-phase 
oscillation unstable, in agreement with the predictions of the BdG analysis. 
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VI. CONCLUSIONS 



In this work, we have performed a detailed study of the statics and dynamics of matter- 
wave dark sohtons in a cigar-shaped Bose-Einstein condensate (BEC) confined in a double- 
well potential. The latter, was assumed to be formed as a combination of a usual harmonic 
trap and a periodic (optical lattice) potential. For our analysis, we have adopted and used 
an effectively ID mean-field model, namely a Gross-Pitaevskii (GP)-like equation with a 
non-cubic nonlinearity, which has been used successfully in other works to describe dark 
solitons in BECs in the dimensionality crossover regime between ID and 3D. 

As a first step in our analysis, we studied the existence and stability of both the one- 
and multiple-dark soliton states. In particular, starting from the respective linear problem, 
we have used the continuation with respect to the chemical potential (number of atoms) 
to reveal the different branches of purely nonlinear solutions, including the ground state, 
as well as excited states, in the form of single or multiple dark solitons. We have shown 
that the presence of the optical lattice is crucial for the existence of nonlinear states that 
do not have a linear counterpart: these are actually asymmetric dark soliton states, with 
solitons located in one well (rather than in both wells) of the double-well potential, and 
do not exist in the non-interacting — linear — limit. We have systematically studied 
the bifurcations of the various branches of solutions and showed how each particular state 
emerges or disappears for certain chemical potential thresholds; the latter were determined 
analytically, in some cases, by using a Galerkin-type approach, in very good agreement with 
the numerical results. For each branch, we have also studied the stability of the pertinent 
solutions via a Bogoliubov-de Gennes (BdG) analysis. This way, we have discussed the role 
of the anomalous (negative energy) modes in the excitation spectra and revealed two types of 
instabilities of dark solitons, corresponding to either purely imaginary or genuinely complex 
eigenfrequencies; the latter were found to occur due to the collision of an anomalous mode 
with a positive energy mode and are associated with the so-called oscillatory instabilities of 
dark solitons. 

We have also adopted a simple analytical model to study small-amplitude oscillations of 
the one- and multiple-dark-soliton states. Particularly, equations of motion for the single and 
multiple solitons were presented in the weakly-interacting limit (where the model becomes 
the usual cubic ID GP equation), and modified them to take into regard the dimensionality of 
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the system (i.e., the effect of transverse directions). The analysis of these equations of motion 
led to the determination of the characteristic soliton frequencies, which were then compared 
to the eigenfrequencies found in the framework of the BdG analysis; the agreement between 
the two was generally found to be very good. Perhaps even more importantly, the simple 
physical (i.e., particle) picture we adopted allowed us to reveal the role of the optical lattice 
strength (i.e., the height of the barrier in the double- well setting) in the stability of the dark 
solitons, as well as its role on the coupling between solitons. Particularly, it was found that 
sufficiently strong barriers lead to instability of symmetric and asymmetric soliton states, 
while it results in an effective decoupling of multi-solitons: the coupling between neighboring 
solitons as described by an effective repulsive potential, becomes negligible for sufficiently 
strong barriers. We should also note that generally the states possessing a dark soliton at 
the center of the trap were found to be most strongly unstable, due to the emergence of a 
progressively larger (the larger the nonlinearity) imaginary eigenfrequency. 

Finally, we have performed systematic numerical simulations, based on direct numerical 
integration of the quasi-one-dimensional model, to investigate the dynamics of dark solitons 
in the double-well setup. We studied the manifestation of instabilities, when present, and 
found in all cases, a very good agreement between the predictions based on the BdG analysis 
and the numerical results. Instabilities have been observed to manifest themselves through 
the exponential divergence of the soliton from its initial location (when associated with 
an imaginary pair of eigenfrequencies) or the oscillatory growth (when associated with a 
complex quartet of eigenfrequencies). 

An interesting direction for a future study would be a similar analysis in a higher- 
dimensional setting: this would refer to vortices confined in double-well potentials (or more 
complex "energy surfaces" such as quadruple-well potentials 
BECs. Such studies are presently in progress. 



71| ) in higher- dimensional 
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